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Abstract 

How does connectivity impact network dynamics? We address this question by pro- 
viding a widely-applicable link between network characteristics on two scales. On the 
global scale we consider the coherence of overall network activity. We show that such 
global coherence can often be predicted from the local structure of the network. To 
characterize this local structure we introduce motif cumulants, which measure the devi- 
ation of pathway counts from those expected in a minimal probabilistic network model. 
Importantly, this link between global dynamics and local architecture is strongly af- 
fected by heterogeneity in the network connectivity. The relationship can be strength- 
ened if the network is divided into connectivity-based subpopulations. Extensions of 
the theory relate higher-order (i.e., beyond pairwise) statistics of network-wide dynam- 
ics and more complex patterns of local connectivity. 



*,''' These authors contributed equally. 
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1 Introduction 



From genetics to neuroscience to the social world, networks of stochastic dynamical systems 
are ubiquitous. The architecture of these networks is complex: irregular but far from random, 



with an overrepresentation of motifs and hubs Bonifazi et al. 



et al. 2011 Milo et al. 2004 Larimer and Strowbridge 2008 . At the same time, networks 



produce complex patterns of collective dynamics Pecora and Carroll 1998 Strogatz 2000 



2009 


Song et al. 


2005, 


Perin 



Rinzel and Ermentrout 1989 . The joint activity of pairs and groups of nodes is frequently 



quantified using correlations and cumulants Renart et al. , 2010 Pernice et al. 2011, 2012 



Lindner et al. 2005 Schneidman et al. 2003 . These measures relate closely to information 



coding and propagation in networks Averbeck et al. , 2006, Salinas and Sejnowski, 2000 



But despite significant progress Pernice et al. , 2011 2012 Renart et al. 2010 Ginzburg and 



Sompolinsky, 1994, Sejnowski, 1976 , relating network activity and connectivity remains a 
challenge. 

Here, our aim is to connect highly local features of a network's connectivity - the frequency 
of small connection pathways, or motifs - with global measures of dynamical coherence. Our 
results hold for any stochastic network where node interactions can be described via linear 
response. These include linearized SDEs (Ornstein-Uhlenbeck) and shot noise processes 



Gardiner , 2009 on networks; our findings for pairwise correlations also extend to integrate- 



Hawkesl p71[ [Pernice et al.[ |2011| ). 



and-fire neurons Trousdale et al. , 2012| and linearly interacting point processes (Hawkes 
models 



We generalize a framework introduced in |Hu et ah 2012, Pernice et al. 2011 to char- 
acterize network structure by the prevalence of certain pathways that connect subsets of 
nodes, quantified by what we call motif moments and motif cumulants. Motif moments 
give a normalized count of pathways of a given length (order). Motif cumulants describe 
deviations in pathway counts at each order from those in a minimal probabilistic model 
of network architecture. Intriguingly, we find that the widely-occuring models of clustered 
networks are highly complex by both measures: high order structure statistics contribute 
significantly to correlated activity. We discuss the origin of this complexity and demonstrate 
a method to resolve it. Finally, we show that the link between global network dynamics and 
local connectivity extends to higher order: n-th order correlations can be predicted from the 
statistics of pathways connecting n nodes. 



Network models Many stochastic networks with linear interactions take the form 

y,{t) = x,{t) + Ai *^Wijyj{t). 



The activity of the i^^ node, yiit), is determined by a (stochastic) baseline activity, Xi{t), 
together with filtered network input. The weight matrix, W, describes the connectivity. To 
focus on network architecture, we make several simplifications: (1) equal-strength connec- 
tions, when present, i.e. W = for an adjacency matrix W°, (2) homogenous nodal 
dynamics, so that Ai{t) = A(t), and (3) statistically uniform (across i) processes Xi(t). These 
assumptions can be relaxed (cf. |Hu et aH|20l2] ). 
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Let y(t) = {yi{t),--- ,yN{t)V and x(t) 
transform of Eq. ([T]) in matrix form a^ 



,xiy(t))'^, and write the Fourier 



(2) 



Here y{uj) = J exp{—2Tciujt){y(t)—Ei [y{t)])dt, x(co') is defined similarly, and A^u) = J^{A){u) 
J^^exp{-2mut)A{t)dt. If the spectral radius '^{A{u)W) < 1, then Eq. implies 
y = (I — y4W)^^x, yielding the following relation between the matrices of cross-spectra: 



Sy{Uj) 



(I-i*W)-^S^(a;)(I-iW' 



T\-l 



(3) 



Denote by a* complex conjugation without transposition. In Eq. ([s]), we used Sj^ 
E [y*y^]; S^.(a;) = E [x*x^], and the Wiener-Khinchin relation Sy{u) = J'(Cy(r)) between 
the cross-spectra and the covariance matrix (Cy(r))jj = E[(yj(t) — Fj[yi{t)]){yj{t + r) — E[yj(t + r)])]. 
Note that under our homogeneity assumptions Sa;(a;) = Sx{uj)1 is a scalar matrix. 

As an illustration, consider Ornstein-Uhlenbeck (OU) processes, which have been used 



widely to model biological networks Pedraza and van Oudenaarden, 2005 Tomioka et al. 



2004, Lestas et al. , 2008, Warren et al. , 2005 . For such processes. 



y = -Ay(t)+Wy(t)+^(t), 



(4) 



where ^{t) is a column vector of white noise processes (here assumed to be independent), 
^I, and T sets the intrinsic timescale of the nodes. Eqs. (2]|3) hold exactly for the 



A 



OU system, with A{uj) = (r ^ + 2Tiiu) ^ and x(w) = A{u)^{uj). An analog of Eq. ^ 



holds approximately in models of neuronal networks [Trousdale et al. 2012 , and exactly for 
Hawkes processes Hawkes, 1971 . Hence our findings about pairwise covariances extend to 
these systems. 



2011 



2012 


Hu et al. , 


2012 


Pernice et al. 



(5) 



n,m=0 



Here, we have normalized by S^iuj) and the unitless quantity on the right measures network 
coherence (we can use this quantity to derive the approximate average correlation coefficient 
- see Supplementary Information (SI)). 

Each term in the sum represents a contribution to the cross spectrum from a differ- 
ent path, or motif, in the network. Consider the second order term AA*(WW'^)ij = 
w^l^P W°;,W°;,. Here w^\A\'^W^,yV'jf, gives the contribution to the spectrum of cells i 
and j of common input (if both connections are present) from node k. In general (W"(W"^)™)jj 
represents the contribution of (n, m) motifs - motifs consisting of two directed chains ema- 
nating from a single apex. Each chain terminates in nodes i and j after traversing n and m 

"'^ Fourier transforms of stochastic processes should be interpreted as taken over finite windows of the 
processes: When computing spectral statistics, scale by 1/T, and take the limit as T — > oo. See Laing and 
LOTdl[2009l 
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nodes, respectively. A (0,m) motif is a chain of length m from node i to node j. The same 
node can appear multiple times in a motif. 

Fig. [l] illustrates this decomposition in a network of two mutually inhibiting nodes. The 
cross-covariance between the nodes, (Cy(r))i2 = J-'~^(Sy(a;))i2, is shown in Fig. [l]^a), while 
contributions of low order motifs to (C2^(r))i2 are given in Fig.[T|^b). As motif order increases, 
corresponding contributions decrease in magnitude, but increase in width. The asymmetry 
of a contribution increases with the asymmetry of the associated motif, i.e. the difference 
between n and m in an {n,m) motif: Compare the contributions of the (1,2) and (0,3) 
motifs. A graphical decomposition of the circuit into the first few (n, m) motifs is shown in 
Fig. Iltc). Since the network is recurrent, the expansion in Eq. ([s]) does not terminate. 




(c) (0,1) (1,2) (0,3) 




Figure 1: (a) The cross-correlation function of two mutually inhibiting nodes, (b) Contri- 
bution of first- and third-order motifs to the cross-correlation function in (a), (c) Graphical 
decomposition of the network showing motifs whose contributions are given in (b). 



Cumulants, moments, and network-wide coherence We now show how averaged 
measures of network coherence depend on network structure. For concreteness - but without 



loss of generality Trousdale et al. , 2012 - we consider integrals of covariance functions over 
long time windows by evaluating all spectral quantities at a; = 0. We indicate this by 
dropping the dependence on u. Moreover, we compute network-wide averages of coherence. 
Denoting by (X) the empirical average of the entries of matrix X, we obtain from Eq. ([s]) 

oo 

n,m=0 

(6) 



N 

n,m=0 



4 



Here the motif moment, = (W''"(W°^)™')/A^"'"'""^ ^, is the empirical probabihty of 

observing an (n,m) motif in the network [Pernice et al. 2011 Hu et al. , 2012 



Note that 



Mn.m = Mm,n- We define /x„n = Mn, and let /Xq n = 1- 



Using Eq. (|3|, we follow Hu et al. , 2012 and obtain an alternate expansion in terms of 
motif cumulants 



Hu et al. , 2012 



Sx 



n+m 



Kir, 



N 



(7) 



The relation between the motif cumulants, K„,m; of Eq. ([T]) and the moments /x„ ^ is similar 

that between cumulants and moments of a random variable. Let C{n) = ■ ■ ■ , nt} : Yli "^i — ^ 

be the set of all compositions {ordered partitions) of n. Then 



{jii,-- ,nt}eC(n) \j=l 

/t-l 

tJ'n,m = E I 11'*^ 

{m,-- ,nt}£C{n) \i=l 
{mi,--- ,ms}£C{m) 




(9) 



In evaluating these terms, we take (n!=i '^nj = 1 if t = 1 and likewise if s = 1. 

The construction of motif moments from cumulants has a familiar interpretation: esti- 
mating the probability of a joint event from the probability of its constituents. Fig. [2|^a) 
demonstrates this for two example motifs. Each term in the diagrammatic expansion arises 
from a cumulant. The first order term gives the probability assuming independence of each 
connection, and the subsequent terms give corrections from excess occurrences of second, and 
then higher-order, submotifs. Hence, each motif cumulants, K„,m captures "pure" higher or- 
der connectivity statistics. For an explicit expression for the motif cumulants Kn,m, see the 
SI and 



Hu et al. 2012 



As an illustrative example, we consider a special case of the two-cluster stochastic block 
model Wang and Wong 1987 Daudin et al. 2008 Litwin-Kumar and Doiron 2012 . Such 



networks are comprised of two subpopulations, each of size iV/2, and associated with a 
constant Sj. The connection probability between nodes in subpopulation i and j is pij = SjS 



J' 



and hence pu = s~ within subpopulation i. While fixing the overall connection probability 
p, the difference between si and $2 describes the degree of clustering in the network. When 
S2 = (si = 2y/p), the only connections are between nodes in the first subpopulation. The 
case Si = $2 = a/p corresponds to an Erdos-Renyi (ER) network. 

Expansions (|6| and ([T]) deliver complementary information about the role of network 
paths in generating coherent activity. Truncating Eq. ^ at order /cmax, so that n+m < fcmax, 
gives the covariance due to paths up to length fcmax- A similar truncation of Eq. ([T]) gives 
the covariance of paths of all orders, with frequencies predicted from paths up to order fcmax- 
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OOOOO O O 0_ kmax 

co" Lo" ID r-" oo" of o" 



Figure 2: (a) The probability of observing two example motifs decomposed into motif cu- 
mulants of the graph, (b) The magnitude of motif cumulants (dashed lines) and moments 
(solid lines) for stochastic block networks: (n, m) motifs with n > m were grouped first 
by n + m and then arranged by increasing n. (c) Approximations of average covariances 
using motif moments (Eq. (|6|, solid lines) and cumulants (Eq. ([t]), dashed lines) up to order 
fcmax- Exact values (direct evaluation of Eq. ([3])) are labeled by crosses. Throughout, we use 
stochastic block models of size 1000 and connection probability 0.2 but with different degree 
of clustering (lighter shades are more clustered). 
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Fig. ^h) shows that motif cumulants decay more rapidly than motif moments. Thus, predic- 
tions of correlations based on cumulants rather than motifs will be more accurate at a given 
order. Second, both motifs and cumulant terms of orders > 1 are larger for more heteroge- 
neous networks. This suggests that heterogeneity in a network's architecture introduces a 
dependence of correlations on longer paths and more complex network structures. 

Fig. ^c) illustrates that the difference is important. At any given order cumulants yield 
more accurate estimates of network coherence. Intriguingly, higher order motif cumulants 
are needed to capture average covariance for other heterogeneous networks. We found that 
networks generated using the Barabasi-Albert (BA) model behaves similarly to a highly 
clustered network Prettejohn et al. , 2011| (see SI), a point we return to below. 



These results agree with our intuition about when local connectivity is a good predictor 
of global dynamics. ER networks are homogeneous, and the most local network statistic - 
connection probability - fully determines graph structure. Thus, for nearly ER networks, 
low order motif cumulants predict dynamical coherence. On the other hand, highly clustered 
networks feature heterogeneous probabilities of linking different nodes. As a consequence, 
the frequency of large motifs depends jointly on multiple connectivity statistics, and cannot 
be obtained accurately from frequencies of lower order motifs exclusively. We next identify 
a novel way of grouping nodes to tame these effects and reestablish the link between local 
connectivity and global coherence. 

Subpopulation cumulants Subsets of nodes can be grouped into classes, or subpopula- 
tions, that share features of dynamics or connectivity. We focus on the latter possibility, and 
characterize each subpopulation by its own motif statistics. Specifically, for b subpopula- 
tions, /x„ „ becomes ab x b matrix of moments. Entry p, q is the empirical probability of an 
(n, m) motif with end nodes belonging to populations p and q, respectively. Motif cumulants 
Kn,m are defined by a recursive relationships similar to Eqs. (8]|9), with products interpreted 



as matrix multiplications (See SI Eqs. (25, 26)). 

We can extend Eq. ([T]) to approximate average covariances using motif statistics within 
and across subpopulations (Eq. (15) in the SI). Fig. [sj^a) illustrates this approach using the 
stochastic block model. If the two sets of nodes used in the generation of the graph define the 
two subpopulations, first order motif cumulants alone predict average correlations perfectly 
(compare with Fig. |2](b)). 

Importantly, the subpopulation approach also works when there is no obvious grouping 
of nodes. As an example we consider finite networks obtained using the BA model. These 
are highly heterogeneous networks with power law degree distributions. If we order the 
nodes according to the sum of in- and out-degree, two subpopulations can be formed from 
nodes with degrees above and below a given threshold value. Fig. shows that result- 
ing predictions of network correlations are significantly improved if the threshold is chosen 
properly. 

With an increasing number of subpopulations, b, we expect the accuracy of predicted 
correlations to increase. As b increases, so does the number of subpopulation statistics at 
our disposal. Thus we have a tradeoff between the resolution in motif statistics at each order 
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and the minimal motif size needed to achieve a given accuracy. As expected, more detailed 
information about local connectivity leads to a better description of global coherence. 




Figure 3: Approximations of average covariances using the subpopulation cumulant approach 
and truncating at order fcmax- Exact values obtained by evaluating Eq. ^ are labeled by 
crosses, (a) stochastic block model networks of Fig. |2] (the colors are the same) divided into 
two subpopulations with differing connectivities; (b) BA network divided into two subpopu- 
lations according to different thresholds on the sum of in- and out-degrees (different colors, 
see inset). 



Higher order correlations We can similarly relate higher order statistics of a network's 
dynamics to its architecture. Our results apply to continuous-valued processes; extensions 
to higher-order coherence for interacting point processes are nontrivial and will be tackled in 
future work. Consider again stationary stochastic processes y{t) = {yi(t), . . . ,?/Ar(t))-^. The 
f^th Qj.(jgj. cross-covariance function is defined using joint cumulants of random variables. 



c;i:p(n,...,^-i) 



f^iyhit),yi2it + n),--- ,?/i^(t + r^.i)) . (lo) 

A generalization of the Wiener-Khinchin theorem relates the Fourier transform of the higher 
order cumulant to the polyspectra S^^*^ "** defined via the Fourier transform of the pro- 
cesses 



Brillinger 



1964 



At third order, we have the bispectrum J^(C*|3j 



s;^[3r(^i,^2) 



E [y*_^{uji + uj2)yi2{(^i)yi3{(^2)] ■ Using Eq. (|2])|^we can generalize Eq. ^ to obtain the bis- 
pectrum for the processes y in terms of that for x and the transfer matrix P = (I — AW)^^: 



s;i3r(^i'^2) 



P:,,,(u;i +u;2)P..,.(u;i)P.3,3(^2)S^f ^'^1 



UJi,UJ2) 



:ii) 



Jl J2 J3 



Expanding P = X]^o(^^)" ^^^"^ leads to an expression for third order polyspectra 

analogous to Eq. (§: 



(Sy[3])/'S'x[3] 



^ oo 



iV2 



:i2) 



n,m,l=0 



^Replacing Gaussian white noise in the OU process with "Poisson kicks", i.e. considering a shot noise 
process, yields non-zero Sa;[3]) 
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(for simplicity, we again set uji = uj2 = 0) ■ Assuming homogeneity, the bispectrum for the un- 
coupled network Sx[3] = Sx[3]lx[3] is a diagonal tensor. More general structure, as well as com- 
mon input, can be treated similarly Hu et al. , 2012 . Here, the motif moments describe con- 



tributions of 3-branch motifs (Fig. 





^' 1,1,1 (hv^^^ l^^f^^,^ ^1^1,1 f^^f^^,^ ^1,1,1 

Figure 4: Cumulant decomposition of a three-branch motif. 



A generalization of Prop. 4.1 in Hu et al. 2012 , allows us to extend Eq. ([7]) to predict 



network-wide third order correlations in terms of motif cumulants: 

(Sy[3]) 




iW] Hr. 

n=l 



Sx[3] 

1 + 3 ^ {NAw)'+^Ki,^ + {NAw)'+'^^^i^i,^,n \ . (13) 



Lm=l Liri,n=l 



Here we use three branch motif cumulants Kn,m,i which have a similar graphical interpretation 
to their second order counterparts (Fig. |4]). The third order motif cumulants may again be 
recursively defined in terms of the motif moments: 




{ni,-- ,ni}eC(n) \i=l 
{mi,--- ,mi,}eC(m) 
{«!,-, /,}GC(0 



(14) 



with the sum being over compositions of n, m, / defined as in Eqs. ( jsP ). 

Our observations for second order statistics generalize directly: Higher order statistics of 
local connectivity predict higher order statistics of network wide dynamics (see SI Fig. [6]). 
First, for both stochastic block models and the BA network, correlations depend significantly 
on long interaction paths though the networks, and on higher-order motifs. Second, at a 
given motif order, motif cumulants provide more accurate predictions of correlations. Third, 
the order of motif statistics needed to approximate correlations again increases with network 
heterogeneity. Finally, the subpopulation approach also generalizes to third order and offers 
similar advantages in predicting correlations from lower-oder motif cumulants (see SI Fig. [?]). 
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Conclusion: For a class of networks with dynamics described by linear response - perhaps 
around an operating point determined by the fully nonlinear dynamics - we have quantified 
the link between local network connectivity and global coherence in activity. In identifying 
the key quantities - motif cumulants - we provide a tool to identify what connection fea- 
tures matter most, and least, in determining collective network dynamics. Moreover, these 
cumulants are both efficient and flexible, in terms of their definition and their ability to pre- 
dict activity based on connection data from small subgroups of nodes. This property could 
provide a way forward in experimental settings - as in studies of networks of genes Alon 
or neurons 



2007 



Song et al. , 2005, Perin et al. , 2011 - in which networks are quantified by 



a limited number of edges are measured simultaneously. 
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Supplementary Information for Local paths to global co- 
herence: cutting networks down to size 

Further examples 

Here we provide details of several computational findings referred to in the main text. Each 
addresses the generality and applicability of our results. First, SI Fig. [5] shows that our 
main results contrasting motif moments and cumulants hold for the Barabasi- Albert network 
model, which has significantly more complex structure than the stochastic block models 
studied in Fig. [2] of the main text. 

Next, SI Figs. |6] and [7] present analogous results for third-order correlations in network 
output. Specifically, SI Fig. [6] shows that these third-order correlations depend significantly 
on the details of the underlying graph structure (i.e., the degree of clustering). Moreover, 
this dependence can be efficiently predicted via motif cumulants. SI Fig. [7] demonstrates 
that the subpopulation approaches continue to enhance the accuracy of our predictions - 
if the populations are correctly defined, levels of triplet correlations can be predicted from 
lower-order motifs. 




CO lo" co" t>-" oo" of o" 



SI Fig. 5: Same as Fig. [2[b,c) of the main text but for the Barabasi-Albert model, (a) The 
magnitude of motif cumulants (dashed lines) and moments (solid lines) for a Barabasi-Albert 
model network. Motifs (n, m), n > m are grouped first by n + m and then arranged by 
increasing n. (b) Approximations of average covariances using motif moments (truncating 
Eq. ([6]), solid lines) and cumulants (truncating Eq. ([T]), dashed lines) up to order /Cmax- Exact 
values (direct evaluation of Eq. (|3|) are labeled by crosses: a Barabasi- Albertnetwork of size 
1000 and connection probability 0.01. 
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SI Fig. 6: Same plots as Fig. [2j^b) of the main text and SI Fig. [5[b) but for average third order 
correlations (Sy[3])/S'a;[3]. (a) Approximations using motif moments (truncating Eq. (12) of 



the main text, solid lines) and cumulants (truncating Eq. (13) of the main text, dashed lines) 
up to order fcmax- Exact values (direct evaluation of Eq. (11) of the main text) are labeled 
by crosses. Graphs were generated from stochastic block models of size 1000 and connection 
probability 0.2 but with different degree of clustering (lighter shades are more clustered), 
(b) Approximations using motif moments (solid lines) and cumulants (dashed lines) up to 
order fcmax for a Barabasi-Albert network of size 1000 and connection probability 0.01. 
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SI Fig. 7: Same plots as Fig. |3j of the main text but for average third order correlations 
{Sy[3]) / Sx[3]. Approximations using the subpopulation cumulant approach by truncating at 
order /cmax, the exact values (direct evaluation of Eq. ([TT]) of the main text) are labeled by 
crosses: (a) stochastic block model networks of Fig. |2] of the main text (the colors are the 
same) divided into two subpopulations with differing connectivities; (b) the Barabasi-Albert 
network of SI Fig. |5] divided into two subpopulations according to different thresholds on the 
sum of in- and out-degrees (different colors, see also the inset, which displays the cutoffs). 



Details of numerical results 

Here we provide a detailed description of the computational examples provided in the main 
text and Supplementary Information. This includes all parameters describing the dynamics 
of nodes and connections, and our methods of generating random networks. 
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In Fig. |T| of the main text, we calculated correlations for an OU system (see Eq. Q of 
the main text) with r = 1, ^ having unit intensity, and 



W 



-0.75 
-0.75 



In plots of approximations of average second and third order covariances, i.e. Fig. ^h), 
[s ^a,b) of the main text, SI Figs. 5](b), 6|a,b), [7|^a,b) and 8|b), the parameters A and w are 
chosen so that NAwp = 0.4. Note that the choice of Sx (resp. S^is] at third order) will not 
affect the normalized quantity {Sy)/Sx (resp. (Sy[3])/5'2.[3]), and can be set to 1. 

The Barabasi- Albert networks in Fig. |3|^b), and SI Figs. [5|^a,b) , [6|(b), [7|b) and[8| ^b) are 
generated by a directed Barabasi- Albert model similar to that in [Prettejohn et al!j 2011| . 
One starts with a "core" of Np nodes, randomly connected with connection probability 0.5. 
After that, N — Np nodes are added to the graph. When adding a new node i + it will form 
exactly Np connections with the existing nodes 1, ■ ■ ■ ,i. Those connections are distributed 
among existing nodes according to probabilities that are proportional to the sum of in- and 
out- degree of each node. The direction of the connection, whether into node i + 1 or out 
of node i + is chosen independently with probability 0.5. The code implementing this 
algorithm is available upon request. 



Explaining the effectiveness of the subpopulation motif approach 



Mathematical underpinnings: In Hu et al. , 2012| , we offered an intuition of the im 



provement in the rate of convergence of the motif cumulant approach (Eq. ([7| )) over the 
motif moment approach (Eq. (|6|). Using the definitions u = (1, ■ ■ ■ , 1)'^ / y/N, H = uu"^, 
= I— H, we compared the spectrum of W° and W°0, which are related to the decay speed 
of ^„ „ and Kn^m respectively. W°0 usually has a much smaller spectral radius since W'^0 
essentially removed the largest isolated eigenvalue Amax- This can be seen from W'^0e = 0, 
where e = (1, ■ ■ ■ , 1)-^ is close to the eigenvector associated with the largest eigenvalue in a 
"homogeneous" network (for example Erdos-Renyi network). 

This reasoning can be extended to treat heterogeneous networks, such as the Barabasi- 
Albert network and stochastic block networks considered here. We begin with the observation 
that the eigenvector for Amax is approximately the (in) degree list d (in our Barabasi-Albert 
networks, the in- and out- degrees are highly correlated). This is also the case for stochastic 
block models. In fact, for any adjacency matrix approximately having the form W = ab^ 
(a, b are column vectors), a is the right eigenvector with eigenvalue b"^a. This indicates that 
the ideas that follow should hold for a general class of networks, including for the Chung 



and Lu model Chung and Lu 2002 



We next apply this to subpopulation motifs for heterogeneous networks. Using the same 
definition of as two paragraphs above will not work in eliminating the remapping largest 
eigenvalue of W*^ to 0. This is because W°0d will not be close to 0, since maps only 
constant vectors to 0. The key in the subpopulation approach is to first redefine the 



matrix with a block diagonal structure, as detailed below SI Eq. (27). Such a will map 
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any vector that is piecewise constant over the indices of each subpopulation to 0. As a 
result, for stochastic block models W°0d = 0, exactly when using the natural division into 
subpopulations. 

For the Barabasi- Albert network, this can only be achieved approximately using only 
two subpopulations and corresponding block diagonal 0. A heuristic argument is that min- 
imizing W°0d, and therefore optimize the accuracy of the subpopulation motif approach, 
is approximately achieved by minimizing ||0d||2. ||0d||2 is determined by the difference 
between d and an approximation of d using piecewise constant vectors associated with the 
subpopulations, which will be different when we choose different cut-offs ranks in the degree 
list. These arguments are confirmed by Fig. [sj^a), where the cut-offs indicated by the analysis 
above reliably lie around the empirically observed optimal value (50) shown in Fig. [sj^b) of 
the main text. 



Extensions to more subpopulations: Following the arguments above, one can refine the 
division of subpopulations in the Barabasi-Albert network to further improve the accuracy 
of the cumulant approach. The Barabasi-Albert network can be viewed as a stochastic block 
model in the limit with large number of subpopulation. This can be confirmed by define a 
that maps d exactly to 0. Indeed, the corresponding cumulant approximations almost 
converge immediately at first order as in stochastic block models (SI Fig. |8](b), compare with 
Fig. |3|^a) of the main text). 



(a) (b) - 10-^ 




degree cut-off k, 



SI Fig. 8: (a) L2 norm of the difference between the degree list (normalized) and the piecewise 
constant vector (see text) given by a certain cut-off. Different shades are 4 realizations of 
Barabasi-Albert networks (with same parameters). The legend is the cut-off value that 
achieves the minimum of difference, (b) Same as Fig. |5|^b) for the 4 realizations of Barabasi- 
Albert networks, using = 1 — dd"^ in Eq. ([T]) (dashed lines) for predicting the weighted 
average -^d'^Syd/ S^- The sohd lines are using the moment approach (Eq. ^ of the main 
text) with a similar replacement of (A) by the weighted average -^d-^Ad. 
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Subpopulation cumulant approach 

Eq. of the main text relates average correlations to motif cumulants for a single pop- 



ulation. Through a similar derivation (see Hu et al. , 2012 for details), we arrive at an 
expression for average correlation in terms of multi-population motif cumulants: 




Y^iNAwr^nE + 5^ {NA 

n=l / \ n,m=l 

-1 



(15) 



oo 



I-2_^{NAw)'^EKr, 

m=l 

where D = diag{l/yiVi, ■ ■ ■ A/VM}, E = ^D'^ = diag{A^i/A^, ■ ■ ■ ,Nb/N}, the matrix 
U is given by U = [ui| ■ ■ ■ |ufc], and Uj = (0, ■ ■ ■ , 0, 1, ■ ■ ■ , 1, ■ ■ ■ , 0)^/ y/Nl is the vector 
where the nonzero entries appear only at indices that match one of the nodes in the given 
subpopulation, normalized to unit L2 norm. Here (8^)^ represents a block-wise average 
over entries corresponding to each subpopulation. We note that by taking the appropriate 
weighted average of these block- wise averages, we arrive again at Eq. ([T]) of the main text. 

Similarly, we can express the third order coherence in terms of motif cumulants for 
networks with multiple subpopulations: 

^ / 00 \ / 00 \ / 00 

(S,[3])/5,[3] = ^ I - Y.(NAwy^iE ® I - Y^iNAwr^mE ® I - 5^(iViu;)"K„E 

\ 1=1 J \ m=l J \ n=l 

00 00 

Y{NAwy+^{Ki,^,.+Ki^.^m + l^.,l,n.)+ Y iNAwy+"'+^Ki^^,^ 
l,m=l l,m,n=l 

(16) 

Here E[3] = diag[3]{A^i/A^, ■ ■ ■ ,Nb/N} is a diagonal tensor. The motif cumulant kz,™,- is a 
tensor defined as 

i'^i^^.r^' = ]^^IT^(DU^K,e,K^UD)„;3/^5. (17) 

This corresponds to a two branch motif, with two endpoints in population a and /5, and 
apex root in population 7. Similarly, Kz,.,m and K.^i^m are simply transpositions of the tensor 
ni,m,-, i-e. (Kz,m,.)"^'^ = {i^i,;mT'^^ = i''^- ,i,mV°'^ ■ Eq. may be verified using Prop, [sj 
which is stated and proven below. 



Note the block average of a N x N x N tensor A may be written a^ (A)b = (DU^'' 



[DU ) (g) (DU ) • A. Applying this averaging identity to Eq. (11) gives that 



(S,[3])b/5x[3] = (DU^)®(DU^)®(DU^)-(P®P®P-I[3]) = (DU^P)®(DU^P)®(DU^P)-I[3]. 



This block average can again be further combined to get the overall average (Eq. ( 13 ) of the 
main text). 



•^Here P (8) Q (g) R • A is understood as the tensor product of matrices P (g) Q (g) R acting on a tensor A, 
that is, (P (g Q (g R • AY^'' = P^QJR^AP'J'-. 
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Explicit expressions for motif cumulants 



Motif cumulants for second order correlations, single population: Here, we will 
prove the following explicit expressions for k„ and Kn^m which are equivalent to the recursive 
definition in Eqs. (8][9) of the main text: 



Kn = y (w°ew° • • • ewo 



N' 

1 

N' 



n factors of W 



-u^ [(wer^wju 



AT" 



u^W^u, 



N' 



y (w°ew° ■ ■ ■ ewo e w°^ew°^ ■ ■ ■ ©w"^^ 

n.+m+l / ^ ^ ^ ^' 



]\fn+r 



n factors of W" 
\n-l 



U^ (We)""^ W0W^ (0W^) 



m factors of W''^ 
j^x »n— 1 



U 



(19) 



JSJn+m 



-u^w^ewlu, 



where 



i)^/yiv, H = uu^, e = I - H. 



and u 



We will assume Eq. (18) is true and demonstrate the proof of Eq. (19). A nearly identical 



(simpler) proof verifies Eq. (18). First, recalling that fi 
straightforward to verify that 



1 



u^(W°)"(W°^)'"u. 



it IS 



(20) 



'"^ ]\[n+m 

Substituting I = + H between every subsequent appearance of the adjacency matrix W° 
gives 



^u^ [W\& + H)] W°(e + H)W°^ [(e + H) W°^] u. (21) 



m— 1 



N' 



By expanding across all sums of + H except the central one (between the terms W'', W°"^), 
and noting that there is an obvious bijection between a pair of compositions of the integers 
n and m, i.e., {rii, . . . ,nt} G C{n), {mi, . . . ,ms} G C(m), and a term of the form 



t-i 



n 



1=1 



[W^^(0 + H)W:,J 



n(Hw^ 
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we may write (using H = uu'^) 
1 



1 



U 



T 



n «H) 



{ni,...,nt}GC{n) Li=l 
_ {mi,...,ms}eC(m) 



[w^je + H)w^^u] 



j=i 



u 



JSJn+n 



E 



{ni,...,nt}eC{?i) 
{mi,...,ms}eC(m) 



t-1 



n (-"w: u) 



i=l 



Kwf.,(e + uu^)w:.u] 



n 



E 



n (^""wf..u 



{jii,.. .,nt}eC(n) L«=l 
{mi,...,ms}eC(m) 



TT ( —n^W u] 

7 = 1 ^ 



(22) 



If t = 1, we define the product [111=1 (W^.H)] = I. 

We now prove Eq. (19) by induction, assuming Eq. (18) holds. First, when n = m = 1, 
the only compositions are trivial (i.e., tti = = {1}). Equating in this case the right-hand 
sides of Eq. ^ of the main text and Eq. ( 22 ) gives that 

2 



Ar2 



Since Eq. ( 18 ) for n = 1 gives that 



we have that Eq. (19) holds for n = m = 1. Next, assume Eq. (19) is true for all (p, q) such 
that p < n and q < m or p < n and q < m. That is, in these cases. 



K, 



- — u^Wju (by Eq. ((18])) and '^p,. - ^ 



Making the corresponding substitutions in Eq. (22), the only term we have not accounted 
for in matching the right-hand side of Eq. (22) to that of Eq. ^ of the main text are the 
terms corresponding to the pair of compositions {n}, {m}. In Eq. ^ of the main text, the 
corresponding terms are 



while in Eq. (22), the terms take the form 

Jn+m 11 m \ n i \ jam 



u'W^(0 + uuMWlu 



(24) 
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where the second equahty follows from the inductive assumption. Equating Eqns. ( 23p4 ) 
gives that 

1 



J\Jn+ra 



-u^w^ewiu, 



which is exactly Eq. (19), completing the inductive proof. 



Motif cumulants for second order correlations, multiple subpopulations: Eq. (15) 



gives a representation of average correlation in terms of motif cumulants under the assump- 
tion of multiple subpopulations. The recursive decomposition relationships for subpopulation 
motif moments and cumulants are 



{rii,...,nt}6C(ra) 



't-1 



,j=l 
't-l 



nt 



{rM,...,nt}eC(ra) \i=l 
{mi,...,ms}eC(m) 



(25) 
(26) 



Here E = d\a.g{Ni/N, ■ ■ ■ , Nb/N} {b is the number of subpopulations, having sizes Ni, where 
Y^Ni = N) is inserted between each motif cumulant matrix multiplication and yields the 
appropriate weighted sums for the interpretation of the terms A^„^ and K„,m as probabilities 
(specifically, scaling by E is multiplication by the probability of selecting cells from respective 
populations at "breaks" in the motifs). Moreover we have an explicit expression for motif 



cumulants similar to Eq. (19): 



Du^ w°ew° ew° e w°^ew°^ ew°^ ud, 



(27) 



n factors of W 



factors of 



where, again, D = diagjl/v^iVi, ■ ■ ■ , l/^/Nh}, = I — H, H = UU"^, the matrix U is given 
by U = [ui| ■ ■ ■ |ufc], and Uj = (0, ■ ■ ■ , 0, 1, ■ ■ ■ , 1, ■ ■ ■ , /y/Wi is the vector where the 
nonzero entries appear only at indices that match one of the nodes in the given subpopulation. 



normalized to unit L2 norm. This formula can be proved identically to Eq. ( 19 ) by considering 
the relationship between j^Y)~^ ^,^^^Y)~^ and ^D~^K„^mD~^ and noting that E = ;^D~^. 



Motif cumulants for third order correlations: The explicit expressions of motif cu- 



mulants in Eqs. (19, 27) also generalize to higher orders. For example, for third-order motif 
cumulants, 



The tensor $ is defined in Prop. [T] 



^E(Wn);(wi)^,(wf)^$^^^ 



(28) 
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For the third order multipopulation case 



fDU^W, 



(29) 



J\fl+m+n-2 ^ ' ' ' ■ 

Here $ is the block diagonal tensor defined in Prop. [3] The summation of a is over all 
subpopulations, and the block tensor $a is defined in Prop. [3j 



Proof of the third order cumulant approach Eq. (13) (main text) 



In Eq. (13) of the main text, we expressed average third-order correlation in terms of motif 



cumulants. Here, we provide the proof of this relationship. We will first prove a proposition 



and Eq. ( 13 ) follows directly. 



Proposition 1. Let H be the rank-1 orthogonal projection matrix generated by the unit 
N -vector u = (1, ■ ■ ■ , 1)/VN, H = uu^, = I - H, tensor ^'^^ = 5'^^ - &^u''/VN - 
&^u^ /^/N — G^''u^/y/N — u^u^u^ /^/N . is the Kronecker delta. We use Einstein sum- 



mation convention and P* 



Ui . For any N x N matrix K, let 
&K, P = (I - K)-^ 



n factors of K 

// the spectral radii \&(K) < 1 and \I^(K©) < 1, then 

oo 



l-5^u^K,u l-5^u^K„u l-^u^K„u 



1=1 

oo 



m=l 



n=l 



l,m=l l,m,n=l 

Proof. Note that, 

= ^/NuiUjUkPlPiPriuPu'^u'' /Vn + BP'^u'/Vn + QP'-uyVN + G^'^mV ViV + $ 



(30) 



pqr\ 



{u,Py){UjPlu,){UkP';Ur) + 3(M,P;0P''P>,)(«fePr 



+VNu,u,UkPlPiP':.^P''- 



p q 



Pu)^ + 3(u' Pu)(u^PeP^u) + y/NuiUjUkPlP{P^r^P'^'' . (31) 

Using Prop. 4.1 i n |Hu et al. [2012 gives alternate forms for the first two terms on the right 
hand side of Eq. (31): 



u'^Pu = ( 1 - 5^ u^KiU 



u^pep^u 



1=1 




(32) 



u^K,eKlu 
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To rewrite the third term, we use a combinatorial lemma similar to that in Hu et al. , 2012 . 



Lemma 2. Let {xi}i>i, {ym}m>i, {z„}n>i, {wimn}i,m,n>i be sequences which converge ab- 
solutely when summed, and also satisfy 



oo 
1=1 



< 1, 



oo 
m=l 



< 1, 



n=l 



< 1. 



Then, 



E E 



i,j,k=l {h,...,lp+i}&C{i) 
{r>ii,...,m,+i}eC{j) 
{ni,...,nr+i}eC(fc) 



oo / oo 

E E^-. 

i=0 \ 1=1 



.s=l 



oo / oo 



A=l 



\u=l 



E E«. 

j=0 \m=l 



E E^ 

k=0 \n=l 

where the sum over {ni, . . . , Ur+i} E C{k) is as defined in the proof of Eqs. (l^TO). 



oo / oo 



oo 

Wlmn 

\Ln,m=l 



Similar to proofs for second order covariances given in Hu et al. , 2012 , using Lemma 2 



and expanding P as a power series in K, we can show that the third term in Eq. (31) may 
be rewritten 



M,M,MfcP;P^^P,^$P'"^= l-^u^K„u J2 w.w,Wfe(K/);(Kj^,(Kj^$^«M . (33) 



n=l 



\Lm,n=l 



Substituting Eqs. ( p2p3| in to Eq. gives Eq. (pOl) after an analytical continuation 
argument similar to that used for the second order analog of Eq. (13) of the main text, 
as shown in |Hu et aL 2012 . 

□ 



Finally, to get Eq. ( 13 ) of the main text, let K = AwW^ in Prop.ll , note that 2UiUjUkSy^^^ 
(Sy[3]), and use the explicit formulas of motif cumulants Eq. (18),([l9|) and (28). 



Proof of the third order subpopulation cumulant approach SI Eq. (16) 



The subpopulation cumulant approach for third order correlation is backed up by the fol- 
lowing proposition. 

Proposition 3. Let U = (ui| ■ • ■ |u;,), and Uj = (0, ■ ■ ■ , 0, 1, ■ ■ ■ , 1, ■ ■ ■ , 0)'^ /y/Wi be the 
unit vector having entries of zero outside of positions corresponding to population i. Also, 
define the 3-dimensional diagonal tensor 

D[3] = diag[3]{l/v^,---,l/v^}, 
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H = UU'^, and = I — H. We refer to the diagonal elements o/Dp] as D^. The tensor $ 
is given by 

b 

a=l 

where 5'^^^ is the Kronecker delta, and (I[3])*-''^ = (5*-^^. Here (0a)*-' := 0a = 0*'' ifi,j are 
both indices in population a or otherwise, is the block of corresponding to population 
a and = Note that $ has a block diagonal structure $ = Ylia^^, where 

is defined as the block of ^: {^^aY^^ = ^^^^ if hj^k are all indices in population a or 
otherwise. Over each block, has a similar expressions as $ in Prop. [1} We use Einstein 
summation convention and A* = Aij = A^^ , = Ui. For any N x N matrix K, let 

Kn = (K0)"^ K = K0K ■ 0K , P = (I - K)^ 

n factors of K 

// the spectral radii \I/(K) < 1 and "^(KB) < 1, then 
(U^P) ® (U^P) ® (U^P) ■ I[3] 

I - J2 U^KzU I ® I I - XI U^K„U j ® I I - J] U^K„,U 

y 1=1 J \ m=l J V n=l y 

oo oo \ 

D[3] + ^i^rn + Yl (U^K,) ® (U^K„) ® (U^K„) ■ $ , (34) 



',m=l Lm,n=l 



where 

[Ri^mT^^ = {lJ^Ki&,KllJ)^^D, + (U^Kz0^K^U)a^D^ + (U^Ki0aK^U)^^D„. (35) 

Proof. The proof parallels exactly the single population case Prop. [Tj Using the definition 
of we have 

(U'^P) ® (U'^P) ® (U^P) • I[3] 
= (U'^PU) ® (U^PU) ® (U^PU) • D[3] + Q 

+(U'^P) (g) (U^P) (g) (U^P) ■ (36) 

where the tensor Q is given by 

Q«/37 = ^(u^P0;,p^u)"^(U^PU)^Da + (U^P0aP^U)"^(U^PU)^Da 

A 

+ (U^P0aP^U)^^(U^PU)"Da. (37) 

Again, the third term is a new type not presented in second order correlation theory. 
It turns out that Lemma |2] holds for this multi-population case where the terms involved 
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are now matrices as opposed to scalars, and scalar multiplications become tensor products. 
Applying the lemma gives, via standard rules for tensor multiplication, 



(U'P) ® (U'P) ® (U'P) ■ $ 

oo \ / oo \ / oo ^ 

1=1 / \ m=l ) \ n=l 

oo 

(U^Ki) ® (U^K^) ® (U'^K„) ■ $ 

\l,m,n=l 

The rest of the proof is identical to the single population case. □ 
Approximating the average correlation coefficient 

The average covariance across the network, (Cy) can be used to approximate the average 



correlation coefficient p^^^ = J2ijPij/-^^y Pij — I \J ^y^v ■ Here we describe two such 
approximations. First, when the uncoupled units have equal mean activities, and the pertur- 
bation from recurrent coupling is weak, the diagonal terms C^' will be close to the uncoupled 
values, Cx- In this case, subtracting diagonal terms from the average, we have that 

P ^\Cx n) N-1 ^^^^ 

In a second, more accurate approximation, we assume permutation symmetry between 
nodes (within a population). Also, in our networks, self-connections are allowed and occur 
with the same probability as other connections. These will lead to identical (marginal) 
distributions for each entry in the covariance matrix Cy, excepting the diagonal entries which 
are shifted by a constant of due to each unit's own uncoupled variance (this corresponds 
to the term proportional to I in expansion Eq. ([s])). This suggests that Cy has the form 



Cy ^ CJ + cInn, 

where Inn is the N x N matrix of ones, and c is a constant. Using the diagonal entries of 
this matrix to normalize, we obtain 

^ Cx + {Cy)-CjN- ^ ^ 

These approximations can also be generalized to networks consisting of multiple sub- 
populations. For example a network with population A and B, we have analogs of Eqns. ( 38p9 ) 
for each of the average covariance within and between two populations A and B), 



avg 

Paa 



[Cy)AA 1 \ 

a naJ'na- 



pTb ^ {Cy)AB/Cx, (40) 
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